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The  numerical  solution  techniques,  including  explicit,  point-implicit,  and  fully  implicit 
schemes,  used  in  a  new  quasi  two-dimensional  procedure  for  the  transient  solution  of  real 
fluid  flows  in  system  lines  and  volumes  are  presented.  The  procedure  is  coupled  with  a  real- 
fluid  properties  database  so  that  both  compressible  and  incompressible  fluids  may  be 
considered  using  the  same  code.  The  procedure  has  been  implemented  in  Matlab/Simulink® 
as  well  as  Fortran95  to  allow  for  application  on  a  wide  variety  of  computer  platforms.  The 
computational  efficiency  of  the  various  numerical  methods  is  discussed  to  aid  in  selection  for 
specific  applications.  Results  for  the  transient  flows  of  gaseous  nitrogen  and  water  in  a 
simple  pipe  network  are  presented  to  demonstrate  the  capability  of  the  current  techniques 
and  the  unsteady  flow  physics  that  can  occur  in  system  lines. 
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Nomenclature 

=  cross-sectional  flow  area 
=  perimeter 

=  relative  acceleration  between  inertial  and  absolute  frames  of  reference  (including  gravity) 
=  Courant-Friedrichs-Lewi  number 
=  speed  of  sound 
=  total  energy  =  p  (e  +  w2/ 2) 

=  fluid  internal  energy 
=  flux  of  mass,  momentum,  and  energy 
=  static  fluid  enthalpy 
=  cell  index 

=  temporal  order  of  accuracy  for  Si  mu  I  ink  ode  15s  solution  methodology 
=  minor  loss  coefficient  =  2Ap/(pif) 

=  index  variables 
=  static  pressure 
=  heat  flux  input 

=  source  vector 
=  time 

=  fluid  state  vector 
=  absolute  velocity 
=  relative  velocity 
=  fluid  cell  volume 
=  distance  along  solution  domain 
=  potential  height  above  sea-level 
=  fluid  density 
=  shear  force  (friction) 
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I.  Introduction 

THE  simulation  of  unsteady  flows  of  non-idealized  fluids  in  system  lines  and  volumes  is  of  interest  for  a  wide 
variety  of  applications  and  industries,  including  air-conditioning  systems,  water/steam/oil  piping  networks, 
refinery  systems,  gas-turbine  secondary  flow-path  and  cooling  networks,  and  liquid  rocket  engine  propellant  lines. 
Several  approaches  for  simulating  the  dynamic  behavior  of  such  fluid-transmission  lines  have  been  reported.  The 
lumped-analysis  approach  treats  a  flow  passage  as  a  series  of  fluid  control  volumes  that  conserve  mass  and  energy 
linked  by  flow  resistance  elements  that  compute  the  flow  between  the  volumes1.  While  this  approach  does  conserve 
momentum  in  a  quasi-steady  sense  at  the  flow  resistances,  the  unsteady  momentum  term  in  the  governing  equations 
is  omitted.  Although  so-called  “continuity”  waves  can  be  captured  using  this  approach,  neglecting  these  terms  leads 
to  an  inability  to  capture  true  “dynamic”  waves  required  for  simulating  such  phenomenon  as  water-hammer  and 
pressure  surge2.  Another  approach  utilizes  the  Method  of  Characteristics3,  which  can  be  applied  to  hyperbolic 
partial  differential  equations.  The  governing  equations  for  fluid  flow  are  compatible  with  this  method  and  it  has 
been  used  for  simulation  of  fluid  transmission  lines4.  While  the  unsteady  momentum  terms  are  retained  using  this 
method,  other  problems,  particularly  at  the  boundaries  of  components,  make  it  difficult  to  apply  to  a  modular 
system-level  simulation  tool.  Modal  methods  have  also  been  used  when  solution  in  the  frequency  domain  is 
possible5.  This  technique  represents  the  pressure  and  velocity  distributions  in  the  flow  domain  as  a  sum  of  an 
infinite  series  of  mode  shapes,  similar  to  a  Fourier  series  solution.  While  this  method  does  present  an  elegant  and 
efficient  method  for  simulating  idealized  flows  (e.g.  incompressible,  inviscid,  laminar,  etc),  the  addition  of  turbulent 
flow,  real-fluid  properties,  heat  transfer  and  phase  change  complicate  the  application  of  the  method  and  reduce  its 
attractiveness. 

Since  unsteady  phenomena  such  as  wave  dynamics  play  an  important  role  in  the  operation  and  testing  of 
systems  that  contain  fluid  lines,  a  method  that  captures  these  transient  effects  is  required.  The  development  of  a 
quasi  two-dimensional,  unsteady,  two-phase  flow  solver  with  heat  transfer  and  real-fluid  properties  using  standard 
finite-difference/control-volume  solution  methods  is  the  subject  of  the  present  effort.  Such  a  solver  is  geared 
towards  modeling  the  dynamic  behavior  of  fluid-filled  lines  and  passages  (i.e.  the  solution  domain  is  much  larger  in 
one  spatial  dimension  than  in  the  others)  accounting  for  the  effects  of  changing  cross-sectional  area.  In  addition,  the 
solver  must  be  suitable  for  use  as  a  module  in  larger  system-level  transient  simulations  of  hydraulic  and  pneumatic 
systems,  so  the  solution  method  must  be  computationally  efficient.  The  following  sections  describe  the  modeling 
approach,  numerical  methodologies  and  test  cases  that  have  been  utilized  during  the  development  of  this  model. 
Results  are  then  shown  for  transient  pipe  flow  of  both  nitrogen  and  water  as  a  demonstration  of  the  numerical 
capability  and  fluid  physics  that  can  be  captured  with  the  current  procedure. 

II.  Approach 

The  model  developed  here  represents  fluid  lines  and  flow  passages  where  the  length  of  the  domain  is  much 
larger  than  the  domain  hydraulic  diameter  so  that  a  quasi  two-dimensional  (2D)  flow  assumption  is  valid.  For  these 
types  of  components,  flow  separations  and  non-axial  velocities  are  minimal,  hence  the  quasi-2D  assumption  is  valid. 
The  solver  is  targeted  to  the  commercial  Simulink®  dynamic  simulation  software  package  from  the  MathWorks  for 
integration  into  a  larger  suite  of  modules  developed  for  simulating  various  systems.  Simulink  '  was  selected  since  it 
offers  a  wide  range  of  capabilities,  over  a  dozen  robust  differential  equation  solvers,  extensive  documentation  and 
technical  support,  a  modem  graphically-based  modeling  paradigm,  an  existing  user  community  across  many 
disciplines,  and  commercially-funded  code  development  and  maintenance.  A  Fortran95  code  using  more  traditional 
solution  methods  is  also  being  developed  in  parallel  to  provide  verification  test  cases  for  the  Simulink  '  module. 

The  solver  is  being  developed  to  account  for  varying  flow  area,  friction,  minor  losses  (e.g.  bends  and  fittings), 
real-fluid  and  two-phase  flow  effects,  gravity  and  acceleration,  heat  transfer,  and  the  capability  to  produce  unsteady 
and  steady-state  solutions.  Fluid  properties  are  obtained  from  the  REFPROP  fluid  property  database6  available  from 
the  National  Institute  for  Standards  and  Testing  (NIST).  This  database  utilizes  state-of-the-art  Equation  of  State 
models  to  fully  describe  fluid  properties  over  a  wide  range  of  thermodynamic  conditions,  including  liquid,  vapor, 
mixed  phase  and  supercritical  fluid  regimes.  Properties  that  completely  define  the  fluid  thermodynamic  state,  as 
well  as  transport  properties,  are  available  as  a  function  of  any  two  thermodynamic  parameters.  Validated  fluid 
models  for  over  80  pure  fluids  and  over  180  fluid  mixtures  are  available  in  the  database.  The  database  is  accessed 
through  a  suite  of  Fortran77  subroutines  that  are  linked  to  the  flow  solver  and  are  used  to  obtain  fluid  equation  of 
state  model  parameters. 

Fluid  friction  and  heat  transfer  are  modeled  as  source  terms  in  the  governing  equations.  This  approach  allows 
the  flow  to  be  modeled  as  one-dimensional  and  facilitates  computational  efficiency.  Friction  (i.e.  viscous  and  minor 
losses)  and  heat  transfer  coefficients  are  obtained  from  suitable  correlations  between  the  flow  variables  and  the 
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source  terms.  The  focus  of  this  paper  is  to  document  the  governing  equations  and  numerical  techniques  used  in  the 
baseline  solution  of  single-phase,  transient  flow  problems  in  the  absence  of  heat  transfer  and  potential  energy. 


III.  Governing  Equations 

The  governing  equations  consist  of  the  quasi  two-dimensional  continuity,  Navier-Stokes  and  energy  equations 
simplified  to  model  viscous  effects  and  general  acceleration  in  a  non-inertial  frame  using  source  terms. 
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where  p  is  the  density,  u  is  the  absolute  velocity  in  the  inertial  frame,  ur  is  the  relative  velocity  in  the  non-inertial 
frame,  a,  is  the  gravitational  and  relative  acceleration  in  the  non-inertial  frame,  E  is  the  total  energy,  p  is  the 
pressure,  H  is  the  stagnation  enthalpy,  V  is  the  local  cell  volume,  A  is  the  local  line  cross-sectional  area,  dx  is  the  cell 
line  length,  Ap  is  the  perimeter,  Tu,  is  the  wall  shear  stress,  and  K  is  the  local  minor  loss  coefficient.  The  wall  shear 
stress  can  be  written  in  terms  of  a  friction  factor,  f  which  is  a  function  of  the  local  Reynolds  number  and  the  wall 
surface  roughness.  For  the  viscous  flow  examples  presented  below,  the  Churchill  correlation7  was  used  to 
determine  single-phase  friction  factors  using  the  equation: 
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where  Re  is  the  local  Reynolds  number  based  on  hydraulic  diameter.  The  minor  loss  coefficient,  K,  may  be 
specified  at  various  locations  within  the  solution  domain  to  model  pressure  loss  at  elbows,  bends,  valves,  and  sudden 
cross  sectional  area  changes. 

Closure  between  the  fluid  pressure,  density,  and  energy  is  performed  with  an  equation  of  state  for  real  fluids 
provided  by  the  NIST  REFPROP  suite  of  thermodynamic  routines.  This  allows  for  the  solution  of  compressible, 
incompressible,  and  two-phase  flows  using  the  same  computational  procedure. 

IV.  Numerical  Methods 

Equations  (1)  and  (2)  may  be  solved  with  a  variety  of  numerical  schemes  that  are  second-order  accurate  in 
space  and  time  for  “steady”  or  transient  flow.  The  solution  schemes  implemented  here  include  implicit,  point- 
implicit,  and  explicit  time-marching  procedures.  Multiple-grid  convergence  acceleration  may  be  used  for  steady- 
flow  problems  as  well  as  during  the  inner  iteration  of  time-accurate,  point-implicit  solutions  to  dramatically  reduce 
the  required  number  of  computational  iterations.  Multiple  numerical  schemes  have  been  developed  using  the 
current  procedure  to  provide  flexibility  and  robustness  when  solving  different  flow  conditions.  All  of  these 
numerical  schemes  use  a  node-centered,  second-order  spatially  accurate,  central-differenced  discretization.  The 
primary  variables  contained  in  the  U  vector  of  Eq.  (2)  are  located  at  the  nodes,  which  are  typically  spaced  equally 
between  the  inlet  and  exit  of  each  line  in  the  overall  domain. 

A.  Implicit  Simulink®  Method 

The  Simulink"  version  of  the  flow  solver  is  implemented  using  an  S-Function  programmed  in  ANSI  C.  An 
S-Function  is  a  specialized  program  for  use  in  Simulink®  models  that  contains  specific  subroutines  required  during 
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model  execution.  The  most  important  of  these  subroutines  defines  the  time  derivatives  of  the  flow  states  (dU/dt)  at 
each  computational  node.  So  for  N  nodes,  this  subroutine  must  define  3N  derivatives,  one  for  each  element  of  the 
state  vector  U  at  each  node.  The  state  fluxes  at  each  cell  center  are  computed  as  discussed  below  (section  IV  B),  and 
the  fluxes  are  distributed  to  the  neighboring  nodes.  Once  the  time  derivatives  of  the  system  states  are  computed, 
standard  Simulink®  solvers  are  used  to  integrate  the  derivatives  and  compute  the  time  history  of  the  model  states. 
Equations  (1)  and  (2)  may  be  solved  iteratively  in  a  sequential  manner.  One  Simulink®  solver  utilized  here,  odel5s, 
treats  the  discretized  form  of  Eqs.  (1)  and  (2)  as  a  sequential  set  of  ordinary  differential  equations  and  solves  them 
iteratively  with  a  Newton  iteration  according  to  the  Numerical  Differentiation  Formula8: 


-At 


f+s 


=  0 


(4) 


where  k  is  the  temporal  order  of  accuracy  of  the  scheme.  In  the  current  investigation,  both  the  temporal  and  spatial 
accuracy  are  second  order.  Another  Simulink®  solver  used  here,  ode23s,  is  a  linearly  implicit  formula  for  stiff 
systems  based  on  a  modified  Rosenbrock  method  that  is  especially  effective  with  crude  tolerances8.  Both  of  the 
Simulink®  solvers  covered  here  are  discussed  in  detail  in  Reference  8. 


B.  Explicit  Method 

Explicit  time-marching  of  Eqs.  (1)  and  (2)  may  also  be  performed  using  a  Lax-Wendroff  control-volume 
scheme,  as  described  by  Ni9.  In  this  procedure,  time-marching  the  primary  variables,  U,  corresponds  to  an  update  in 
time  at  each  node  according  to  a  second  order  temporally  accurate  Taylor  series  formula: 
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The  Ni  scheme  is  a  finite  volume  integration  method.  The  scheme  is  applied  in  a  three-step  process.  First,  the 
change,  At  rJU/dt,  (the  second  term  on  the  right-hand-side  of  Eq.  (5))  at  the  center  of  each  computational  cell  is 
approximated  as: 
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The  time-step  size  is  determined  by  the  Courant-Friedrichs-Lewi  condition  i.e. 
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where  CFL  is  the  stability  number  with  a  typical  value  of  0.7  and  c  is  the  local  speed  of  sound. 

The  correction,  0.5  At2  52U/3t2,  (the  third  term  on  the  right-hand-side  of  Eq.  (5))  to  the  points  i  and  i+1  are 
determined  from  time  derivatives  of  Eq.  (6).  These  corrections  are  then  added  or  subtracted  to  the  cell  center 
change  given  by  Eq.  (6)  to  obtain  the  time -rate  contributions  to  the  adjacent  nodes  from  each  cell: 
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The  primary  variables  at  the  nodes  are  then  updated  by 

&jt=@ut  +^Ui)cellt  +  (dUM  +d2UM)ce!lj 


(11) 


U"+1  =  U”  +  SUi  (12) 

The  CFL  stability  number  that  sets  the  time-step  size  is  set  to  0.7  due  to  the  explicit  nature  of  this  scheme.  For 
steady  flows  or  the  inner  iteration  of  the  point-implicit  technique  described  below,  a  multiple-grid  acceleration 
procedure  described  by  NiQ  is  used  to  greatly  reduce  the  number  of  time-steps  to  reach  convergence.  This  allows  the 
computational  efficiency  of  this  scheme  to  be  competitive  with  the  implicit  procedure  described  above. 

C.  Point  Implicit  Method 

A  point-implicit  procedure  can  be  made  up  from  either  the  implicit  or  explicit  schemes  described  above.  For  this 
procedure,  the  time-rate  changes  in  Eqs.  (4)  or  (5)  are  assumed  to  be  in  pseudo-time.  A  true  time  derivative  of  the 
primary  variables  is  added  to  the  right-hand  side  of  Eq.  (1).  This  true  time  derivative  is  derived  from  stored 
solutions  of  the  primary  variables  at  k+1  previous  time-steps,  where  again,  k  is  the  temporal  order  of  accuracy 
desired.  An  inner-iteration  is  then  used  to  drive  the  right-hand-side  to  zero.  The  advantage  of  this  approach  is  that 
steady-flow  acceleration  techniques  may  be  used  during  the  inner  iteration  such  as  a  multiple-grid  scheme.  Also, 
this  technique  provides  flexibility  in  performing  multiple  lines  of  a  network  in  parallel.  This  procedure  is  known  as 
the  dual-time-step  technique10  and  has  many  similarities  to  the  Newton  iteration  used  in  the  implicit  procedure. 

D.  Boundary  Conditions 

Boundary  conditions  are  applied  each  time-step  at  the  inlet  and  exit  of  each  line  in  the  network.  These  boundary 
conditions  are  specified  depending  on  the  network  connectivity  and  consist  of  either  a  global  inlet,  global  exit,  or 
junction  boundary  condition.  At  the  global  inlet(s)  to  the  network,  the  stagnation  enthalpy  and  entropy  are  held 
constant  at  levels  determined  from  specified  stagnation  pressure  and  stagnation  temperature  and  the  predicted 
velocity.  The  resulting  density  and  internal  energy  are  determined  from  REFPROP  thus  allowing  for  the 
determination  of  the  total  energy,  E.  At  the  global  exit(s)  to  the  network,  the  specified  exit  static  pressure  is  held 
constant.  The  internal  energy  and  predicted  velocity  along  with  this  specified  pressure  are  used  with  REFPROP  to 
determine  the  density  and  total  energy.  The  prescribed  inlet  stagnation  pressure,  stagnation  temperature,  and  exit 
static  pressure  may  be  specified  as  functions  of  time.  At  line  junctions,  the  contributions  of  the  time -rate  changes  in 
the  primary  variables  contained  in  U  of  Eq.  (1)  are  summed  to  give  the  total  time-rate  change  at  the  junction  node. 
For  each  junction  of  a  given  line,  the  neighboring  line  numbers  and  boundaries  (i.e.  inlet  or  exit)  are  stored  to  make 
the  implementation  of  this  boundary  condition  straightforward. 

V,  Results 

Transient  pipe  flow  simulations  of  gaseous  nitrogen  and  liquid  water  have  been  performed  to  verify  the  accuracy 
and  demonstrate  the  capability  of  the  present  procedure.  A  co-planar  (z  =  0)  two-pipe  system,  shown  in  Fig.  1, 
consisting  of  two  0.127  m  (5  in)  pipes  with  0.0254  m  (1  in)  diameters  was  used  for  these  simulations.  The  surface 
roughness  of  the  pipe,  necessary  for  viscous  flow  simulations,  was  2.54  x  10"6  m  (0.0001  in).  The  inlet  stagnation 
pressure  as  well  as  the  initial  exit  static  pressure  was  held  at  345  kPa  (50  psia).  Thus,  the  initial  velocity  of  the  flow 
in  the  pipe  was  zero.  The  inlet  stagnation  temperature  was  held  at  300  K  (540  R).  At  time  0.001  seconds,  the  exit 
pressure  was  reduced  to  33 1  kPa  (48  psia).  This  sudden  reduction  in  pressure  initiates  a  series  of  expansion  and 
compression  waves  through  the  pipe  system  and  an  increase  in  the  velocity  of  the  flow.  For  inviscid-flow  cases,  the 
pressure  decays  in  the  pipe  to  a  uniform  value  and  the  velocity  grows  asymptotically  to  a  uniform  value.  For 
viscous-flow  cases,  the  pressure  in  the  pipe  decays  but  to  a  non-uniform  distribution  with  the  inlet  pressure  higher 
than  the  exit  due  to  friction  and  minor  losses.  The  velocity  of  the  pipe  flow,  however,  will  once  again  grow 
asymptotically  in  time  until  it  reaches  a  near-uniform  value.  All  simulations  were  run  in  time  until  both  the  pressure 
and  velocity  distributions  through  the  pipe  remained  time-independent.  A  total  of  33  grid  points  with  uniform 
spacing  were  used  along  each  pipe  in  order  to  discretize  the  governing  equations. 
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Figure  1.  Two-pipe  system. 
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A.  Transient,  Viscous  Nitrogen  Pipe  System  Flow 

The  two-pipe  flow  simulation  with  nitrogen  is  meant  to  test  the  current  procedure’s  capabilities  for  modeling 
compressible  flow.  The  explicit  and  point-implicit  numerical  procedures  were  run  in  order  to  verify  the  solution 
accuracy  and  to  demonstrate  features  of  these  solution  techniques.  Figure  2  shows  the  predicted  pressure  and 
velocity  as  a  function  of  time  using  the  explicit  numerical  technique  for  viscous  flow  at  five  locations  along  each 
pipe.  In  this  scheme,  the  time  step  size  is  determined  by  the  minimum  time  step  of  all  the  cells  in  the  overall 
domain.  The  sinusoidal  variation  of  pressure  with  time  is  evidence  of  the  expansion  and  compression  waves 
propagating  through  the  pipe.  The  time-variation  of  the  velocity  corresponds  to  these  waves.  The  pressure  decays 
in  time  until  it  reaches  equilibrium  with  an  overall  pressure  drop  across  the  system  of  5.47  kPa  (0.793  psia).  The 
maximum  pressure  drop  occurs  at  the  exit  of  line-1  and  the  inlet  of  line-2  where  minor  losses  occur  due  to  the  flow 
turning  in  the  90  degree  elbow.  The  velocity  increases  until  it  reaches  a  value  of  65.59  m/s  (215.2  ft/s).  The  friction 
factor  predicted  by  Eq.  (3)  was  0.00188.  The  Reynolds  number  of  the  flow  based  on  line  diameter  was  predicted  to 
be  3.66xl05.  Similar  results  were  obtained  using  the  time-accurate  implicit  Simnlink  solvers  odel5s  and  ode23s. 


Line  1  Line  2 


Figure  2.  Nitrogen  transient  viscous  pipe  flow  using  explicit  numerical  technique. 
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Line  1 


Line  2 


Figure  2.  (cont)  Nitrogen  transient  viscous  pipe  flow  using  explicit  numerical  technique. 


The  solution  upon  reaching  steady-state  may  be  compared  with  an  analytical  solution  to  the  energy  equation. 
For  incompressible  flow,  the  control-volume  energy  equation  may  be  written  as: 
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where  hfriction  is  the  viscous  friction  head,  g  is  the  acceleration  due  to  gravity,  and  z  is  the  potential  height  above  a 
reference  location.  At  the  inlet  to  the  configuration  shown  in  Fig.  1,  the  pressure  is  equal  to  the  stagnation  pressure 
and  the  velocity  is  zero.  At  the  exit  to  line -2,  the  pressure  is  equal  to  the  exit  static  pressure.  The  friction  head  may 
be  written  in  terms  of  the  friction  and  minor  loss  coefficients  as: 

f  \ 


1  friction 

where  L  is  the  pipe  length,  and  Dh  is  the  hydraulic  diameter.  For  the  problem  described  in  Fig.  1  using  nitrogen  as 
the  fluid,  the  analytical  solution  of  Eq.  (13)  results  in  a  steady-state  velocity  of  65.67  m/s  (215.4  ft/s)  and  a  pressure 
drop  of  5.45  kPa  (0.791  psia).  The  predicted  velocity  corresponds  to  a  low  enough  Mach  number  (0.15)  that  the 
nitrogen  gas  can  be  assumed  essentially  incompressible.  The  agreement  between  the  numerical  prediction  and  the 
steady-state  analytical  solution  is  excellent  with  the  minor  differences  being  due  to  the  low  level  of  compressibility 
in  the  gas.  The  numerical  solution  shown  in  Fig.  2  took  5,850  iterations  to  reach  steady-state  over  0.0397  seconds  in 
real  time. 

As  a  check  of  the  transient  response  of  the  computational  solution,  the  period  of  pressure  oscillations  can  be 
compared  to  the  period  expected  for  frictionless  flow  (e.g.  organ  pipe  modes).  The  pressure  waves  inside  the  pipe 
will  travel  at  slightly  less  than  the  speed  of  sound,  where  the  difference  is  a  result  of  friction  losses.  Thus,  the 
pressure  oscillations  should  occur  close  to  the  natural  frequency  of  the  pipe  (c/2L),  or  696Hz  for  this  case.  The 
observed  oscillation  period  is  about  0.00157  seconds,  which  is  slightly  higher  than  the  period  corresponding  to  the 
natural  frequency  of  the  pipe  (0.00144  seconds).  Thus,  the  method  produces  transient  results  in  keeping  with 
physical  expectations. 

Figure  3  shows  the  predicted  results  for  the  same  two-pipe  system  with  nitrogen  using  the  point-implicit 
numerical  procedure.  In  this  procedure,  the  time-step  size  can  be  prescribed  by  the  user  to  be  any  value  that  results 
in  a  desired  solution  time  and  resolution  of  temporal  features  in  the  flow.  For  the  present  case,  the  time-step  size 
was  set  at  0.01  seconds,  which  was  approximately  1,700  times  larger  than  that  required  for  in  the  explicit  procedure. 
The  capability  to  use  a  large  time-step  size  in  the  point-implicit  procedure  allows  for  very  fast  solution  time. 
However,  as  shown  by  comparing  Figs.  2  and  3,  resolution  of  the  high-frequency  waves  and  their  effect  on  the 
pressure  and  velocity  is  lost  as  a  result  of  the  increased  time-step  size.  The  availability  of  different  solution 
algorithms  to  the  overall  procedure  allows  for  rapid  solutions  when  fine  temporal  resolution  is  not  needed  as  well  as 
detailed  solutions  when  high  temporal  resolution  of  waves  is  required. 
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Figure  3.  Nitrogen  transient  viscous  pipe  flow  using  point-implicit  numerical  technique. 


B.  Transient,  Viscous  Water  Pipe  System  Flow 

A  similar  set  of  simulations  were  performed  using  the  present  procedure  with  water  to  demonstrate  the  capability 
to  predict  transient  flows  of  incompressible  fluids.  The  main  critical  difference  with  incompressible  fluids  is  that  the 
speed  of  sound  and  wave  speeds  are  much  larger  than  for  compressible  fluids.  This  makes  the  use  of  explicit 
numerical  techniques  quite  time-consuming.  Resolution  of  the  detailed  unsteady  sinusoidal  transients  shown  above 
with  the  explicit  numerical  technique  is  possible,  but  requires  large  compute  times.  The  effects  of  the  high- 
frequency  waves  are  not  always  of  interest,  so  that  numerical  techniques  that  can  take  large  time-steps  and  reduce 
solution  time  are  often  attractive. 

Figure  4  shows  the  transient,  point-implicit,  viscous  solution  of  water  through  the  previously  described 
configuration  (Fig.  1).  Note  that  the  time  to  reach  steady-state  for  water  is  much  longer  than  that  required  for 
nitrogen  due  to  the  incompressibility  of  the  fluid.  The  size  of  the  time  step  used  in  the  point-implicit  solution  was 
0.1  seconds.  The  time  step  size  required  in  the  corresponding  explicit  solution  procedure  was  approximately  2xl0'6 
seconds  or  50,000  times  smaller  than  that  used  in  the  point-implicit  solution.  As  a  result,  approximately  1.5  trillion 
time  steps  would  be  required  to  solve  the  3000  seconds  of  time  required  to  reach  steady-state.  Even  with  a  time  step 
size  of  0.1  seconds,  the  point-implicit  scheme  still  resolves  much  of  the  pressure  oscillations  in  the  flow.  The 
predicted  pressure  drop  over  the  two-pipe  system  converges  to  a  steady-state  value  of  5.59  kPa  (0.811  psia)  which 
again  is  in  excellent  agreement  with  the  analytical  energy  equation  analysis  that  gave  5.61  kPa  (0.814  psia).  The 
velocity  converges  to  a  uniform  steady-state  value  of  4.05  m/s  (13.29  ft/s).  The  predicted  friction  factor  was 
0.00228  and  the  Reynolds  number  based  on  diameter  was  1.06x10s.  The  analytical  energy  equation  analysis 
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resulted  in  a  velocity  of  4.05  m/s  (13.28  ft/s)  demonstrating  that  the  numerical  techniques  are  verified  against 
analytical  solutions. 

Line  1  Line  2 


Time  (s)  Time  (s) 

Pressure 


Time  (s)  Time  (s) 

Velocity 

Figure  4.  Water  transient  viscous  pipe  flow  using  point-implicit  numerical  technique. 

Figure  5  shows  the  predicted  transient  pressure  over  the  time  period  between  0.001  and  0.003  seconds  for  both 
lines  resulting  from  the  Simulink8  implicit  numerical  technique.  Also  shown  in  this  figure  is  the  predicted  transient 
pressure  over  the  same  time  period  predicted  with  the  point-implicit  numerical  scheme.  The  time-step  size 
determined  automatically  in  the  Simulink"  implicit  technique  (i.e.  2xl0'5  seconds)  was  approximately  10  times 
greater  than  that  which  would  have  been  used  in  the  explicit  technique.  The  same  time-step  size  of  2xl0"5  seconds 
was  used  in  the  point-implicit  scheme.  As  with  the  point-implicit  procedure,  the  fully  implicit  technique  has  the 
advantage  of  being  able  to  use  very  large  time  step  sizes.  This  capability  affords  fast  turnaround  of  solutions  with 
the  possibility  of  some  loss  in  temporal  resolution.  However,  both  of  the  implicit  and  point-implicit  schemes  also 
allow  for  small  time  step  sizes  similar  in  magnitude  to  those  used  by  the  explicit  numerical  technique  in  order  to 
resolve  the  transient  details  and  frequencies  in  the  flows  if  desired.  This  ability  to  use  a  wide  range  of  time  step 
sizes  allows  for  maximum  flexibility  in  obtaining  the  desired  solutions. 

Comparison  of  the  unsteady  pressure  between  the  implicit  and  point-implicit  numerical  techniques  in  Fig.  5 
shows  that  the  predicted  frequency  of  the  unsteady  waves  is  essentially  the  same.  The  observed  pressure  oscillation 
period  of  about  0.000341  seconds  compares  well  with  the  expected  period  of  0.000338  seconds  (2,957Hz)  based  on 
frictionless  flow  analysis.  The  predicted  amplitudes  of  the  unsteady  waves  of  the  Simulink8  implicit  procedure  are 
greater  than  those  predicted  by  the  point-implicit  scheme.  In  addition,  the  Simulink8  implicit  procedure  predicts 
more  harmonics  in  the  unsteady  solution  than  predicted  by  the  point-implicit  scheme.  These  differences  are  likely 
due  to  the  dispersion  error  resulting  from  the  adaptive  time-step  used  in  the  Simulink"  implicit  procedure. 
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Figure  5.  Transient  pressure  comparison  of  numerical  techniques  for  water  transient  viscous  pipe  flow. 


VI.  Conclusion 

A  new  baseline  procedure  has  been  developed  for  the  numerical  solution  of  transient  quasi  two-dimensional  flow 
in  system  lines,  networks,  and  volumes.  This  new  procedure  has  been  implemented  in  both  Matlab/Simulink"  and 
Fortran95.  A  variety  of  numerical  solution  algorithms  are  available  to  allow  for  various  needs  depending  on  the 
desired  solution  time  and  frequency  resolution.  Currently,  explicit,  point-implicit,  and  implicit  numerical  techniques 
have  been  implemented  and  verified.  The  implicit  and  point-implict  techniques  are  useful  for  obtaining  rapid  turn¬ 
around  of  solutions  at  the  expense  of  some  loss  in  temporal  resolution.  The  explicit  technique  is  useful  for 
predicting  the  high-frequency  content  of  transient  flow  fields.  Transient,  viscous  solutions  using  these  numerical 
techniques  have  been  presented  and  verified  for  a  two-pipe  system  with  a  closed-form  analytical  solution  to  the 
energy  equation. 

Agreement  between  the  numerical  and  analytical  solutions  based  on  energy  balance  analysis  and  pressure 
oscillation  frequency  are  excellent  demonstrating  that  the  present  procedure  for  the  transient  solution  of 
compressible  or  incompressible  fluids  in  the  absence  of  heat  transfer  is  valid.  Unsteady  pressure  and  velocity 
through  the  two  lines  in  the  system  are  shown  as  a  demonstration  of  the  transient  physics  predicted  by  the  procedure. 

Work  on  the  quasi  two-dimensional  flow  solver  codes  is  continuing.  Future  efforts  will  focus  on  implementing 
two-phase  flow  and  wall  heat  transfer  capabilities,  modularization  for  use  in  higher-level  system  simulations,  and  an 
approach  for  parallel  solution  on  multi-processor  configurations. 
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